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Abstract. In this work we present numerical study of a trapped Bose-Einstein condensate 
perturbed by an alternating potential. The relevant physical situation has been recently realized 
in experiment, where the trapped condensate of ®^Rb, being strongly perturbed, exhibits the 
set of spatial structures. Firstly, regular vortices are detected. Further, increasing either the 
excitation amplitude or modulation time results in the transition to quantum vortex turbulence, 
followed by a granular state. Numerical simulation of the nonequilibrium Bose-condensed 
system is based on the solution of the time-dependent 3D nonlinear Schrodinger equation within 
the static and dynamical algorithms. The damped gradient step and time split-step Fourier 
transform methods are employed. We demonstrate that computer simulations qualitatively 
reproduce the experimental picture, and describe well the main experimental observables. 


1. Introduction 

As is well known, trapped Bose atoms in equilibrium at zero temperature form Bose-Einstein 
condensate (BEC) in the ground state. This system has been thoroughly studied theoretically 
as well as experimentally. For details we refer to the review articles and books m El El m 
and references therein. If the condensate is slightly moved from equilibrium it usually exhibits 
collective excitations. Such dynamics is typical not only for BEC, but for the variety of quantum 
systems mum- 

Nonequilibrium BEC can be created by injecting into the system energy, which can be realized 
by applying an external modulation field. Two main ways of moving trapped BEC far from 
equilibrium are usually considered. The first one is a resonant excitation of coherent modes 
mm- In this case the applied external perturbation is weak, but the modulation frequency of 
the alternating field is in resonance with one of the transition frequencies between topological 
modes. This method is a bit difficult for experimental realization, requiring fine tuning to 
resonance. The other way is to employ a strong external alternating field superimposed onto the 
stationary trap [inj, where either the modulation time or excitation amplitude determine the 
energy injected into the system. In both these cases, BEC transfers from its ground to excited 
states, demonstrating a set of dynamical regimes. The second approach has been recently 
realized in the experiment mm- 

The experimental realization of nonequilibrium condensates stimulates interest to its 
numerical modeling and computer simulation m- Dealing with a trapped Bose-Einstein 
condensate, the nonlinear Schroedinger equation (NLSE), which, in BEC physics sometimes 


is called the Gross-Pitaevskii equation [T], is the main tool for the theoretical investigation. 
Nowadays, different numerical schemes [laii! have been developed and successfully applied for 
the solution of second order partial differential equations, like NLSE. The methods differ in their 
stability, accuracy, computational cost, and ability of keeping the NLSE original properties. In 
general, all such methods could be classihed as explicit, implicit and finite Eourier transform 
or pseudospectral methods. Each numerical scheme has its own merits and disadvantages, 
determining the accuracy of the NLSE solution. The presence of either random forces or 
potentials, being a prerequisit of nonequilibrium systems, complicates the numerical procedure, 
since the relevant NLSE solution is not smooth. 

In our recent paper m we presented the results of numerical simulation of a trapped 
condensate, under external perturbation, showing the appearance of a granular state, analogous 
to that observed in the experiment. The simulations were based on the solution of the NLSE, 
with the parameters corresponding to the experimental setup. In the present paper we give the 
details of the accomplished numerical modeling. Also, we demonstrate that the chosen numerical 
methods qualitatively describe the typical physical picture, observed in the experiment. 

The paper is organized as follows. In Sec. 1, we briefly review the experiment on the strong 
perturbation of a trapped BEG by external modulating field. Details of computer simulation 
are given in Sec. 2. In Sec. 3, we present the numerically simulated time-amplitude diagram. 

2. Experiment on perturbation of 

In the experiment mm, the ®^Rb atoms (atomic mass m = 1.444 g, scattering length 
tts = 0.557 X 10“® cm), being cooled down to the temperature much lower than = 276 
iiK, form the ground state Bose-Einstein condensate. The total number of atoms in the BEG 
fraction is ~ 1.5 x 10^. Stationary trapping potential has cylindrical symmetry, with the trap 
frequencies: = 27r x 210 Hz, = 27r x 23 Hz. After the stationary BEG state, obeying the 

Thomas-Eermi shape, is achieved, the trapping potential is modulated by an additional external 
alternating field, oscillating with the frequency Hq = 27r x 200 Hz. As a result, the total trap 
potential Vext can be approximated as [12]: 


^ext 




+ [ycos0i — a: sin 01 — < 12(1 — cosHo^)]^ + 
-|-[ 2 :cos 02 + xsin 02 — ^ 3(1 — cos Hot)]^), 


( 1 ) 


where A = 0* = Aj(l — cos Hot), Ai = ^, A 2 = The trapping frequency Ur 

determines the oscillator length that defines the modulation parameters of the potential: 
(di, 62 , S 3 ) = a(2,5, 3) fim /Ir, where a is an effective amplitude of pumping. 

If the external perturbation imposes a certain angular momentum onto the system, like in the 
case of the rotating m or laser-stirred condensate m, the formed vortices are aligned along 
the rotation axis. The number of vortices depends on the modulation frequency. The aligned 
vortices can form a completely ordered lattice. 

But in our case, as is seen from the expression the modulation field used in the 
experiment does not impose a certain rotation axis. This feature makes the main difference 
of our experiment, as compared to the case of rotating BEG. Contrary to the rotated BEG, the 
external perturbation, we use, just shakes the condensate cloud, injecting the energy into the 
trap and thus generating the diverse coherent topological modes [lOj . 

Mathematically, coherent topological modes are stationary solutions of NLSE, including the 
ground state and other excited states. Different modes correspond to rather different spatial 
distributions of the condensate density, forming a variety of spatial shapes. As a result, different 
structures can be generated in a trapped condensate [Eli. 


First modes, generated in the experiment, are quantum vortices. Since the external held 
does not inject a rotation moment, the vortices appear as vortex-antivortex pairs with the 
winding numbers plus/minus one These modes are energetically the most stable. The number 
of created vortices depends on the injected energy, hence can be controlled by either the pumping 
amplitude or by the modulation time. At a hxed amplitude, almost a linear dependence of the 
vortex number from time is observed, until the distance between the vortex centers is less than 
four coherence lengths, when the interaction energy of two vortices is much smaller than the 
vortex energy. When the number of vortices becomes sufficiently large, the transition to random 
vortex turbulence occurs. 

The continuous energy pumping into the BEC cloud creates more vortices than the trap can 
host, as a result of which the vortices start strongly interacting and colliding, which destroys 
the regime of vortex turbulence. Then the system passes to the next dynamical regime, where 
the condensate forms high-density grains surrounded by sufficiently rarified gas. It can be 
demonstrated m that these grains are coherent objects. 

3. Numerical simulation 

Experimental characteristics are found to be in good agreement with analytical estimates [iniiisi- 
Now we aim at accomplishing numerical calculations for the same setup as in the experiment. 
In this section, we focus on the techniques of numerical modeling as applied to the solution of 
NLSE describing the condensate dynamics driven by external perturbations. The corresponding 
numerical simulation is realized in two steps: first, we calculate the parameters of the initial 
conditions describing the stationary state of the BEC in the trap, and then we consider time 
propagation. The latter is performed via the solution of the time-dependent NLSE for the desired 
time interval. In both these cases, the full 3D geometry is considered, and the wave function 
is determined on the 3D space grid. Below, the main points of the accomplished numerical 
simulations are illustrated. 

3.1. Initial conditions 

First, the wave function, describing the condensate stationary state, has to be computed. To 
this end, we consider the eigenproblem: 

^ — Eyi^Pyii ( 2 ) 

where the index n enumerates quantum states and the corresponding wave functions (pn and 
energies In equation ([2|), the nonlinear Hamiltonian H[p) is typical for a trapped BEC: 

+ gN\p\^ , (3) 

Im 

where g = is the interaction strength depending on the scattering length and atomic 

mass m, Vext is the trapping potential. The parameters of the Hamiltonian ([2D exactly 
correspond to the experimental setup. The BEC of ®^Rb atoms with the related Og and m 
are considered. The total number of condensed particles \s N = 1.5 x 10®. The external field 
Vext is represented by the total trap potential in form ([ID. Since we interested in stationary 
solutions, the external perturbation is switched off, so that Vext = Vext{^,t = 0). 

The eigenproblem ([2D can be solved by a variety of numerical schemes, but iterative methods 
are the most convenient. We used the damped gradient step method [Ulllg], where the initial 
guess (jtn is improved by the iteration steps. Mathematically, this can be represented in the form: 

>)<^^] , 


( 4 ) 



where 2? is a damping operator, the symbol O labels the orthonormalization of the set of 
functions (pn at each iteration step marked as i. The Schmidt orthonormalization is usually 
employed. Hamiltonian ([3]) is recalculated after the each iteration step. The main problem 
of the numerical scheme 0 is the proper choice of the damping operator D providing the 
convergence of the iteration procedure. The damped gradient step method operates with the 
damping operator in the form 


V = 


5 

T + Eo ’ 


(5) 


where 5 is a parameter controlling the iteration step size, T = —is the operator of kinetic 
energy. Since kinetic energy can happen to be close to zero, parameter Eq stabilizes the iteration 
process. 

If the procedure converges, the final set of pn, obtained after a number of iteration steps, 
corresponds to the eigenfunctions with the spectra E^- 

After the eigenproblem ([2]) is solved, the total wave function 'k(r) can be written as a 
superposition of the states 


and then normalized to the total number of condensed atoms 


N = 


/ 


|T(r)pdr . 


(7) 


At t = 0, there is no external perturbation, thus the condensate is in equilibrium state 
corresponding to the lowest energy Eq. To take this into account, we apply the condition 


Cq » Cn>0 


( 8 ) 


to the notation ([5|) . 

Representation ([6]) of the wave function corresponds to its expansion over the basis of coherent 
topological modes [9], which allows one to treat non-ground-state condensates [8]. 

3.2. Time evolution 

The static solution is followed by the time propagation of the condensate wave function p{r). 
For this purpose, we solve numerically the nonlinear 3D time-dependent NLSE 

^ [-^V^ + Kxt(r)+5|^(r,t)p]T(r,t) , (9) 

using the obtained 'l'(r) as the initial condition. Similarly to the previous case, the parameters 
of NLSE exactly correspond to the experimental setup. To solve NLSE, the time split-step 
Fourier transform (TSSFT) method [201121j is chosen. The advantage of this numerical scheme 
m is that it is unconditionally stable and keeps the main properties of NLSE, such as mass 
conservation, time reversibility, i.e., time-inverse invariance. Also, the method is second-order 
accurate in both time and space. 

The main idea of the TSSET is in splitting the time-propagation between the kinetic and 
potential components of the total energy. Both these components can be approximated by 
exponents related to the solutions of the linear Schrodinger equation (see Ref. [21] for details). 
Mathematically, the mechanism of time propagation has the following form: 


T(r, f -|- At) ^ exp {iTAt/2) exp {—iVAt) exp (2rAt/2)T(r, t) 


( 10 ) 



where T = —and V = Vext + are the operators of kinetic and potential energy, 

respectively. 

The relation between the exponential operators obeys the Baker-Cambell-Hausdorf theorem, 
which states that is valid if C = A + B + B] + .... Since the TSSFT method 

skips the commutation relations of the kinetic and potential operators, numerical error grows 
with time. The symmetrical form of (|10ll . with respect to the operators, allows to remove the 
hrst commutator term, increasing the accuracy of time propagation [22j . 

The exponential derivative operator is easily calculated in the momentum space rather than 
in the coordinate space, because of the well known relation 

exp = exp(-A:2)T'[T(q,t)] , (11) 

where symbol J- means the Fourier transform. A: is a variable conjugated to q. 

Taking into account expressions (fTUI) and m, the full time propagation of the condensate 
wave function for the time step At is performed in three stages: 

(i) The Fourier transform of the initial function 'I'(r,t) is followed by the half time-step At/2, 
with the operator of kinetic energy in the momentum space. The first stage is completed 
by the inverse Fourier transform that results in the new function 'I'(r,t) ^ T(r,t). 

(ii) Recalculation of the potential operator V with the function T(r,t) and further full-time 

step propagation of 'k(r, t), with the recalculated term V. As a result, the transformation 
T(r,t) —> is obtained. 

(iii) The hnal stage equals the first step, but with the function T(r, t). This completes the time 

propagation, since ^(r,t) + At) . 

The time step At should be small enough to keep the accuracy of time propagation. The 
increment At fulhlls the condition [23] 


< TlFtlT > 


27r 

<< ^ 

At 


( 12 ) 


3.3. Other numerical aspects 

Dealing with a grid-based technique, the parameters of the computation box are extremely 
important. The size of the computational box strictly depends on the trapping frequencies, 
total number of condensed particles and the kind of boundary conditions. The proper size of 
the 3D grid avoids the unphysical truncation of the operated wave function and its reflection 
from the grid boundaries during the perturbation. The later is important, if reflecting boundary 
conditions are employed. In the present simulation, we used the reflecting boundary conditions 
for both, the static and dynamical simulations. In the experiment, the applied external 
perturbation heats up the BEC leading to the transition from condensed to thermal phase for a 
part of atoms. This process is almost regular in time. Consequently, the number of condensed 
particles in the trap permanently decreases. To follow the experiment, we introduced into the 
system a small dissipation replacing the term ih-^^{r,t) in the NLSE (Hj) by (i — 'y)h-^'i>{r,t), 
where 7 is a phenomenological term. The same approach is also used by the other authors [23]. In 
simulations, during the 60 ms of perturbation, approximately 15 percents of the initial condensed 
particles are lost, which is in good agreement with the experiment. Of course, introducing the 
dissipation term 7 does not conserve the initial norm of the wave function ([7|), which implies 
the diminishing number of condensed atoms. 


4. Results and discussion 

Computer simulations reproduce well the experimental results. The applied external 
pertnrbation, modelled by ([I]) creates, first, dipole oscillations of the condensate cloud, followed 
by vortex structures corresponding to vortex-antivortex pairs. The vortex nature of the observed 
structures is confirmed by quantized circulation having the values plus or minus one. The 
circulation is calculated by the method proposed in Ref [25] . Increasing either the amplitude of 
pumping or the excitation time produces more and more vortices, creating the random vortex 
tangle, representing, according to Feynman [2B], quantum turbulence. The next regime observed 
in the simulations is grain turbulence, where continuous trap modulation destroys all the vortices, 
forming high-density grains surrounded by rarihed gas. The typical physical picture, obtained 
in numerical simulations, is represented by the time-amplitude diagram shown in Fig. 1. The 
experimental diagram, published in Ref m, is in qualitative agreement with that obtained in 
the numerical simulation. The area inside the rectangle reproduces the experiment in the best 
way. 



Figure 1. The time-amplitude diagram, as found in numerical simulations. The area inside the 
rectangle is in the best agreement with the experiment. 

Numerical simulations give very good quantitative values for the main characteristics observed 
in the experiment and also found from theoretical arguments. Thus, the critical number of 
vortices Ncr ~ 25, required for creating the regime of quantum turbulence, is the same in 
the experiment, numerical simulations, as well as in theoretical estimates. The grains can 
be considered as coherent objects, since their typical linear size practically coincides with the 
coherence length, as is found in the experiment and in the simulation. 

In numerical simulations, we also observe one more regime [15], arising after the granular 
state, the regime of wave turbulence, when all BEC is destroyed, and turbulence is formed by 
small-amplitude waves. Here we compare the experimental and numerical results, because of 
which the regime of wave turbulence is not included in the diagram, since it has not yet been 
reached in experiment. 
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